Re: [R] creating a color gradient in geom_ribbon

[Ulrik Stervbo]

I haven't tested it  but the first thing I'd look at is scale_fill_gradient.

HTH


Ulrik

Jim Lemon  schrieb am Do., 11. Mai 2017, 07:22:

> Hi Kristi,
> It can be done, but it is messy:
>
> pl = data.frame(Time = 0:10, menle = rnorm(11))
> pl$menlelb = pl$menle -1
> pl$menleub = pl$menle +1
> rg<-0.95
> blue<-1
> plot(pl$Time,pl$menlelb,ylim=range(c(pl$menlelb,pl$menleub)),type="l",
>  lwd=7,col=rgb(rg,rg,blue))
> lines(pl$Time,pl$menlelb,lwd=7,col=rgb(rg,rg,blue))
> rg<-seq(0.9,0.3,length.out=9)
> offset<-seq(0.88,0.08,by=-0.1)
> for(i in 1:9) {
>  lines(pl$Time,pl$menle+offset[i],lwd=7,col=rgb(rg[i],rg[i],blue))
>  lines(pl$Time,pl$menle-offset[i],lwd=7,col=rgb(rg[i],rg[i],blue))
> }
> lines(pl$Time,pl$menle,lwd=6,col=rgb(0,0,blue))
>
> For the ggplot solution, this might work:
>
> ggplot(pl, aes(Time)) +
>   geom_line(aes(y=menle+1), colour=rgb(0.95,0.95,1), width=7) +
>   geom_line(aes(y=menle-1), colour=rgb(0.95,0.95,1), width=7) +
>   geom_line(aes(y=menle+0.88), colour=rgb(0.9,0.9,1), width=7) +
>   geom_line(aes(y=menle-0.88), colour=rgb(0.9,0.9,1), width=7) +
>   geom_line(aes(y=menle+0.78), colour=rgb(0.825,0.825,1), width=7) +
>   geom_line(aes(y=menle-0.78), colour=rgb(0.825,0.825,1), width=7) +
>   geom_line(aes(y=menle+68), colour=rgb(0.75,0.75,1), width=7) +
>   geom_line(aes(y=menle-68), colour=rgb(0.75,0.75,1), width=7) +
>   geom_line(aes(y=menle+0.58), colour=rgb(0.675,0.675,1), width=7) +
>   geom_line(aes(y=menle-0.58), colour=rgb(0.675,0.675,1), width=7) +
>   geom_line(aes(y=menle+0.48), colour=rgb(0.6,0.6,1), width=7) +
>   geom_line(aes(y=menle-0.48), colour=rgb(0.6,0.6,1), width=7) +
>   geom_line(aes(y=menle+0.38), colour=rgb(0.525,0.525,1), width=7) +
>   geom_line(aes(y=menle-0.38), colour=rgb(0.525,0.525,1), width=7) +
>   geom_line(aes(y=menle+0.28), colour=rgb(0.45,0.45,1), width=7) +
>   geom_line(aes(y=menle-0.28), colour=rgb(0.45,0.45,1), width=7) +
>   geom_line(aes(y=menle+0.18), colour=rgb(0.375,0.375,1), width=7) +
>   geom_line(aes(y=menle-0.18), colour=rgb(0.375,0.375,1), width=7) +
>   geom_line(aes(y=menle+0.08), colour=rgb(0.3,0.3,1), width=7) +
>   geom_line(aes(y=menle-0.08), colour=rgb(0.3,0.3,1), width=7) +
>   geom_line(aes(y=menle), colour="blue") )
>
> but I can't test it.
>
> Jim
>
> On Thu, May 11, 2017 at 6:05 AM, Kristi Glover
>  wrote:
> > Hi R Users,
> >
> > I was trying to create a figure with geom_ribbon. There is a function
> "fill", but I want to make the shaded area with a gradient (increasing dark
> color towards a central line, inserted of having a color). Is there any
> possibility?
> >
> >
> > In the given example, I want the colour with "blue" but in a gradient
> (dark=central, light= as goes higher or lower)
> >
> >
> > pl = data.frame(Time = 0:10, menle = rnorm(11))
> >
> > pl$menlelb = pl$menle -1
> >
> > pl$menleub = pl$menle +1
> >
> > ggplot(pl, aes(Time)) +
> >
> >   geom_line(aes(y=menle), colour="blue") +
> >
> >   geom_ribbon(aes(ymin=menlelb, ymax=menleub), fill="blue")
> >
> >
> > Thanks
> >
> > [[alternative HTML version deleted]]
> >
> > __
> > R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
>
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] creating a color gradient in geom_ribbon

[Jim Lemon]

Hi Kristi,
It can be done, but it is messy:

pl = data.frame(Time = 0:10, menle = rnorm(11))
pl$menlelb = pl$menle -1
pl$menleub = pl$menle +1
rg<-0.95
blue<-1
plot(pl$Time,pl$menlelb,ylim=range(c(pl$menlelb,pl$menleub)),type="l",
 lwd=7,col=rgb(rg,rg,blue))
lines(pl$Time,pl$menlelb,lwd=7,col=rgb(rg,rg,blue))
rg<-seq(0.9,0.3,length.out=9)
offset<-seq(0.88,0.08,by=-0.1)
for(i in 1:9) {
 lines(pl$Time,pl$menle+offset[i],lwd=7,col=rgb(rg[i],rg[i],blue))
 lines(pl$Time,pl$menle-offset[i],lwd=7,col=rgb(rg[i],rg[i],blue))
}
lines(pl$Time,pl$menle,lwd=6,col=rgb(0,0,blue))

For the ggplot solution, this might work:

ggplot(pl, aes(Time)) +
  geom_line(aes(y=menle+1), colour=rgb(0.95,0.95,1), width=7) +
  geom_line(aes(y=menle-1), colour=rgb(0.95,0.95,1), width=7) +
  geom_line(aes(y=menle+0.88), colour=rgb(0.9,0.9,1), width=7) +
  geom_line(aes(y=menle-0.88), colour=rgb(0.9,0.9,1), width=7) +
  geom_line(aes(y=menle+0.78), colour=rgb(0.825,0.825,1), width=7) +
  geom_line(aes(y=menle-0.78), colour=rgb(0.825,0.825,1), width=7) +
  geom_line(aes(y=menle+68), colour=rgb(0.75,0.75,1), width=7) +
  geom_line(aes(y=menle-68), colour=rgb(0.75,0.75,1), width=7) +
  geom_line(aes(y=menle+0.58), colour=rgb(0.675,0.675,1), width=7) +
  geom_line(aes(y=menle-0.58), colour=rgb(0.675,0.675,1), width=7) +
  geom_line(aes(y=menle+0.48), colour=rgb(0.6,0.6,1), width=7) +
  geom_line(aes(y=menle-0.48), colour=rgb(0.6,0.6,1), width=7) +
  geom_line(aes(y=menle+0.38), colour=rgb(0.525,0.525,1), width=7) +
  geom_line(aes(y=menle-0.38), colour=rgb(0.525,0.525,1), width=7) +
  geom_line(aes(y=menle+0.28), colour=rgb(0.45,0.45,1), width=7) +
  geom_line(aes(y=menle-0.28), colour=rgb(0.45,0.45,1), width=7) +
  geom_line(aes(y=menle+0.18), colour=rgb(0.375,0.375,1), width=7) +
  geom_line(aes(y=menle-0.18), colour=rgb(0.375,0.375,1), width=7) +
  geom_line(aes(y=menle+0.08), colour=rgb(0.3,0.3,1), width=7) +
  geom_line(aes(y=menle-0.08), colour=rgb(0.3,0.3,1), width=7) +
  geom_line(aes(y=menle), colour="blue") )

but I can't test it.

Jim

On Thu, May 11, 2017 at 6:05 AM, Kristi Glover
 wrote:
> Hi R Users,
>
> I was trying to create a figure with geom_ribbon. There is a function "fill", 
> but I want to make the shaded area with a gradient (increasing dark color 
> towards a central line, inserted of having a color). Is there any possibility?
>
>
> In the given example, I want the colour with "blue" but in a gradient 
> (dark=central, light= as goes higher or lower)
>
>
> pl = data.frame(Time = 0:10, menle = rnorm(11))
>
> pl$menlelb = pl$menle -1
>
> pl$menleub = pl$menle +1
>
> ggplot(pl, aes(Time)) +
>
>   geom_line(aes(y=menle), colour="blue") +
>
>   geom_ribbon(aes(ymin=menlelb, ymax=menleub), fill="blue")
>
>
> Thanks
>
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Solve system of non linear equations using nasted loops

[Bert Gunter]

I haven't gone through your code carefully, but I believe this can be done
in a tiny fraction of the time you are taking by eschewing loops. See
?expand.grid to get started. In general,nested loops should be avoided if
possible.

I  also suggest you spend some time with a good R tutorial to learn how to
write good R code instead of recapitulating C (or whatever) in R. All
languages have their own coding paradigms(R's is functional) and you will
do better in the long run if you spend the time to learn R's.

Cheers,
Bert



On May 10, 2017 1:20 PM, "Santiago Burone" 
wrote:

Hello,

I'm new at R and I would like to use it in order to solve a system of non
linear equations. I have the code that works but im not able to save the
results.



My system has three equations and i would like to solve this using three
nested loops, becouse i need solutions for all the values combinations.





The code im using is this:



library(nleqslv)
#x3 es gamma
#x2 es alpha
#x1 es beta
MA <- c(5, 43600, 4, 38800, 37600, 34400, 31600, 27200, 24400,
2)
MI <- c(1, 21800, 2, 19400, 18800, 17200, 15800, 13600, 12200,
1)
ME <- c(3, 32700, 3, 29100, 28200, 25800, 23700, 20400, 18300,
15000)
DE <- c(0.384900179, 0.19245009, 0.19245009, 0.19245009, 0.19245009,
0.19245009, 0.19245009, 0.19245009, 0.19245009, 0.19245009)
for (i in 1:9) {
for (j in 1:9){
for (k in 1:9){
df.names <- paste("SociedadB",1:10,sep="")

fun <- function(x) {
f <- numeric(length(x)) # read as:
f[1] <- 1*x[3] - (log(ME[1+i]/ME[1])+(1*x[1]+1*
x[2])*(2)*(log((MA[1]-ME[1])/(MA[i+1]-ME[i+1]/(log(0.
19245009/0.384900179))
f[2] <- 1*x[3] - (log(MI[1+j]/MI[1])-(1*x[2])*(
2)*(log((MA[1+j]-MI[1+j])/(MA[1]-MI[1]/(log(0.19245009/0.384900179))
f[3] <- 1*x[3] - (log(MA[1+k]/MA[1])-(1*x[1])*(
2)*(log((MI[1+k]-MA[1+k])/(MI[1]-MA[1]/(log(0.19245009/0.384900179))
f
}
startx <- c(1,1,1)
answers<-as.data.frame(nleqslv(startx,fun))
d.frame <- (answers$x)
assign(df.names[k+20], d.frame)
despejes<- data.frame(SociedadB2, SociedadB3, SociedadB4, SociedadB5,
SociedadB6, SociedadB7, SociedadB8, SociedadB9, SociedadB10)

}
}
}

despejes<- data.frame(SociedadB2, SociedadB3, SociedadB4, SociedadB5,
SociedadB6, SociedadB7, SociedadB8, SociedadB9, SociedadB10)



If im not wrong i should have 1000 different solutions but the data frame
only save 10 results. I tried copying the last part of the code between the
last two brackets (changing the name of the data frame despejes for
despejes 1, 2 and 3 but i only get three data frames with exactly the same
results. I dont know how could i save all the results of every combination.
I mean, i need to solve for i=1, j=1 and the 9 diferent values of k and
save those results, then i=1, j=2 and the 9 values of k..



I know this might be a easy question but i have not been able to find a
solution in the last days.



Thanks in advance if you can help me.


Santiago.
[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Solve system of non linear equations using nasted loops

[David Winsemius]

> On May 10, 2017, at 11:20 AM, Santiago Burone  
> wrote:
> 
> Hello, 
> 
> I'm new at R and I would like to use it in order to solve a system of non 
> linear equations. I have the code that works but im not able to save the 
> results. 
> 
> 
> 
> My system has three equations and i would like to solve this using three 
> nested loops, becouse i need solutions for all the values combinations.
> 
> 
> 
> 
> 
> The code im using is this:
> 
> 
> 
> library(nleqslv)
> #x3 es gamma
> #x2 es alpha
> #x1 es beta
> MA <- c(5, 43600, 4, 38800, 37600, 34400, 31600, 27200, 24400, 2)
> MI <- c(1, 21800, 2, 19400, 18800, 17200, 15800, 13600, 12200, 1)
> ME <- c(3, 32700, 3, 29100, 28200, 25800, 23700, 20400, 18300, 15000)
> DE <- c(0.384900179, 0.19245009, 0.19245009, 0.19245009, 0.19245009, 
> 0.19245009, 0.19245009, 0.19245009, 0.19245009, 0.19245009)
> for (i in 1:9) {
> for (j in 1:9){
> for (k in 1:9){

I got an error on the first pass through the loop.


> df.names <- paste("SociedadB",1:10,sep="")

# seems kind of excessive to repeat that process 9*9*9 times
# Maybe this should be done outside the loop?

> 
> fun <- function(x) { 
> f <- numeric(length(x)) # read as:
> f[1] <- 1*x[3] - 
> (log(ME[1+i]/ME[1])+(1*x[1]+1*x[2])*(2)*(log((MA[1]-ME[1])/(MA[i+1]-ME[i+1]/(log(0.19245009/0.384900179))
> f[2] <- 1*x[3] - 
> (log(MI[1+j]/MI[1])-(1*x[2])*(2)*(log((MA[1+j]-MI[1+j])/(MA[1]-MI[1]/(log(0.19245009/0.384900179))
>  
> f[3] <- 1*x[3] - 
> (log(MA[1+k]/MA[1])-(1*x[1])*(2)*(log((MI[1+k]-MA[1+k])/(MI[1]-MA[1]/(log(0.19245009/0.384900179))
>  
> f 
> } 

# So you need 9*9*9 different functions? 


> startx <- c(1,1,1) 
> answers<-as.data.frame(nleqslv(startx,fun)) 

> d.frame <- (answers$x)
> assign(df.names[k+20], d.frame)

# I don't think that df.names[k+20] would be anything other than NA since it 
only has 10 elements

> despejes<- data.frame(SociedadB2, SociedadB3, SociedadB4, SociedadB5, 
> SociedadB6, SociedadB7, SociedadB8, SociedadB9, SociedadB10)
> 
> }
> }
> }
> 
> despejes<- data.frame(SociedadB2, SociedadB3, SociedadB4, SociedadB5, 
> SociedadB6, SociedadB7, SociedadB8, SociedadB9, SociedadB10)
> 

The error message is:

Error in data.frame(SociedadB2, SociedadB3, SociedadB4, SociedadB5, SociedadB6, 
 : 
  object 'SociedadB2' not found

That's because column names are not first class data-objects in R.


> 
> 
> If im not wrong i should have 1000 different solutions but the data frame 
> only save 10 results. I tried copying the last part of the code between the 
> last two brackets (changing the name of the data frame despejes for despejes 
> 1, 2 and 3 but i only get three data frames with exactly the same results. I 
> dont know how could i save all the results of every combination. I mean, i 
> need to solve for i=1, j=1 and the 9 diferent values of k and save those 
> results, then i=1, j=2 and the 9 values of k..
> 

If you need to save separate results for each combination of i,j,k , then you 
need at a minimum an array (possibly with a further 4th dimension if you are 
storing numerical vectors)  or possibly a list if these items are more 
complicated than an atomic vector.
> 
> 
> I know this might be a easy question but i have not been able to find a 
> solution in the last days. 
> 
> 
> 
> Thanks in advance if you can help me.
> 
> 
> Santiago. 
>   [[alternative HTML version deleted]]

Rhelp is a plain text list although you don't seem to have run into trouble on 
that account. You should however determine how your mail client can be setup to 
just send pain text to this mailing list.


> 
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

David Winsemius
Alameda, CA, USA

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] Solve system of non linear equations using nasted loops

[Santiago Burone]

Hello, 

I'm new at R and I would like to use it in order to solve a system of non 
linear equations. I have the code that works but im not able to save the 
results. 

 

My system has three equations and i would like to solve this using three nested 
loops, becouse i need solutions for all the values combinations.

 

 

The code im using is this:

 

library(nleqslv)
#x3 es gamma
#x2 es alpha
#x1 es beta
MA <- c(5, 43600, 4, 38800, 37600, 34400, 31600, 27200, 24400, 2)
MI <- c(1, 21800, 2, 19400, 18800, 17200, 15800, 13600, 12200, 1)
ME <- c(3, 32700, 3, 29100, 28200, 25800, 23700, 20400, 18300, 15000)
DE <- c(0.384900179, 0.19245009, 0.19245009, 0.19245009, 0.19245009, 
0.19245009, 0.19245009, 0.19245009, 0.19245009, 0.19245009)
for (i in 1:9) {
for (j in 1:9){
for (k in 1:9){
df.names <- paste("SociedadB",1:10,sep="")

fun <- function(x) { 
f <- numeric(length(x)) # read as:
f[1] <- 1*x[3] - 
(log(ME[1+i]/ME[1])+(1*x[1]+1*x[2])*(2)*(log((MA[1]-ME[1])/(MA[i+1]-ME[i+1]/(log(0.19245009/0.384900179))
f[2] <- 1*x[3] - 
(log(MI[1+j]/MI[1])-(1*x[2])*(2)*(log((MA[1+j]-MI[1+j])/(MA[1]-MI[1]/(log(0.19245009/0.384900179))
 
f[3] <- 1*x[3] - 
(log(MA[1+k]/MA[1])-(1*x[1])*(2)*(log((MI[1+k]-MA[1+k])/(MI[1]-MA[1]/(log(0.19245009/0.384900179))
 
f 
} 
startx <- c(1,1,1) 
answers<-as.data.frame(nleqslv(startx,fun)) 
d.frame <- (answers$x)
assign(df.names[k+20], d.frame)
despejes<- data.frame(SociedadB2, SociedadB3, SociedadB4, SociedadB5, 
SociedadB6, SociedadB7, SociedadB8, SociedadB9, SociedadB10)

}
}
}

despejes<- data.frame(SociedadB2, SociedadB3, SociedadB4, SociedadB5, 
SociedadB6, SociedadB7, SociedadB8, SociedadB9, SociedadB10)

 

If im not wrong i should have 1000 different solutions but the data frame only 
save 10 results. I tried copying the last part of the code between the last two 
brackets (changing the name of the data frame despejes for despejes 1, 2 and 3 
but i only get three data frames with exactly the same results. I dont know how 
could i save all the results of every combination. I mean, i need to solve for 
i=1, j=1 and the 9 diferent values of k and save those results, then i=1, j=2 
and the 9 values of k..

 

I know this might be a easy question but i have not been able to find a 
solution in the last days. 

 

Thanks in advance if you can help me.


Santiago. 
[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] creating a color gradient in geom_ribbon

[Kristi Glover]

Hi R Users,

I was trying to create a figure with geom_ribbon. There is a function "fill", 
but I want to make the shaded area with a gradient (increasing dark color 
towards a central line, inserted of having a color). Is there any possibility?


In the given example, I want the colour with "blue" but in a gradient 
(dark=central, light= as goes higher or lower)


pl = data.frame(Time = 0:10, menle = rnorm(11))

pl$menlelb = pl$menle -1

pl$menleub = pl$menle +1

ggplot(pl, aes(Time)) +

  geom_line(aes(y=menle), colour="blue") +

  geom_ribbon(aes(ymin=menlelb, ymax=menleub), fill="blue")


Thanks

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Generating samples from truncated multivariate Student-t distribution

[David Winsemius]

> On May 10, 2017, at 11:02 AM, Czarek Kowalski  wrote:
> 
> Previously I had used another language to make calculations based on
> theory. I have calculated using R and I have received another results.
> My theoretical calculation does not take into account the full
> covariance matrix (only 6 elements from diagonal). Code based on
> theory:
> 
> df = 4;   #degrees of freedom
> sigmas = c(1, 4, 2, 5, 3, 6) # roots of diagonal elements of covariance matrix
> meann = c(55, 40, 50, 35, 45, 30)
> alfa1 = 20; # lower truncation
> beta1 = 60; # upper truncation
> a = (alfa1 - meann)/sigmas;
> b = (beta1 - meann)/sigmas;
> E = meann + sigmas * ((gamma(df - 1)/2)*((df + a^2)^(-(df-1)/2) - (df
> + b^2)^(-(df-1)/2))*df^(df/2))/(2*(pt(b,df)-pt(a,df))*gamma(df/2)*gamma(1/2))

This looks wrong:

(gamma(df - 1)/2)

According to your "theory" (which you have not yet supported with references) 
in the attached pdf file, that should be:

gamma( (df - 1)/2 )

I'm not a statistician and checking your "theoretical" expression for the mean 
of a truncated PDF is not really on-topic for R help. When I look at 
http://www.tonyohagan.co.uk/academic/pdf/trunc_multi_t.PDF I see a lot of 
normalization factors that seem quite different than yours. Perhaps you should 
post any further difficulties to stats.stackexchange.com ?

You may also want to consult:

https://www.jstatsoft.org/article/view/v016c02


Authors:Saralees Nadarajah, Samuel Kotz
Title:  R Programs for Truncated Distributions
Abstract:   Truncated distributions arise naturally in many practical 
situations. In this note, we provide programs for computing six quantities of 
interest (probability density function, mean, variance, cumulative distribution 
function, quantile function and random numbers) for any truncated distribution: 
whether it is left truncated, right truncated or doubly truncated. The programs 
are written in R: a freely downloadable statistical software.

Best;
David.


> E
> 
> 
> Kind regards
> Czarek
> 
> On 9 May 2017 at 23:50, David Winsemius  wrote:
>> 
>>> On May 9, 2017, at 2:33 PM, David Winsemius  wrote:
>>> 
>>> 
 On May 9, 2017, at 2:05 PM, Czarek Kowalski  wrote:
 
 I have already posted that in attachement - pdf file.
>>> 
>>> I see that now. I failed to scroll to the 3rd page.
>>> 
 I am posting
 plain text here:
 
> library(tmvtnorm)
> meann = c(55, 40, 50, 35, 45, 30)
> covv = matrix(c(  1, 1, 0, 2, -1, -1,
  1, 16, -6, -6, -2, 12,
  0, -6, 4, 2, -2, -5,
  2, -6, 2, 25, 0, -17,
 -1, -2, -2, 0, 9, -5,
 -1, 12, -5, -17, -5, 36), 6, 6)
 df = 4
 lower = c(20, 20, 20, 20, 20, 20)
 upper = c(60, 60, 60, 60, 60, 60)
 X1 <- rtmvt(n=10, meann, covv, df, lower, upper)
 
 
> sum(X1[,1]) / 10
 [1] 54.98258
 sum(X1[,2]) / 10
 [1] 40.36153
 sum(X1[,3]) / 10
 [1] 49.83571
 sum(X1[,4]) / 10
 [1] 34.69571  # "4th element of mean vector"
 sum(X1[,5]) / 10
 [1] 44.81081
 sum(X1[,6]) / 10
 [1] 31.10834
 
 And corresponding results received using equation (3) from pdf file:
 [54.97,
 40,
 49.95,
 35.31, #  "4th element of mean vector"
 44.94,
 31.32]
 
>>> 
>>> I get similar results for the output from your code,
>>> 
>>> My 100-fold run of your calculations were:
>>> 
>>> meansBig <- replicate(100, {Xbig <- rtmvt(n=10, meann, covv, df, lower, 
>>> upper)
>>> colMeans(Xbig)} )
>>> 
>>> describe(meansBig[4,])  # describe is from Hmisc package
>>> 
>>> meansBig[4, ]
>>>  n  missing distinct Info Mean  Gmd  .05  .10  
>>> .25
>>>1000  1001 34.7  0.0195434.6834.68
>>> 34.69
>>>.50  .75  .90  .95
>>>  34.7034.7234.7234.73
>>> 
>>> lowest : 34.65222 34.66675 34.66703 34.66875 34.67566
>>> highest: 34.72939 34.73012 34.73051 34.73742 34.74441
>>> 
>>> 
>>> So agree, 35.31 is outside the plausible range of an RV formed with that 
>>> package, but I don't have any code relating to your calculations from 
>>> theory.
>> 
>> Further investigation:
>> 
>> covDiag <- covv*( row(covv)==col(covv) )  # just the diagonal means
>> 
>> Repeat with all zero covariances:
>> 
>>> meansDiag <- replicate(100, {Xbig <- rtmvt(n=10, meann, covDiag, df, 
>>> lower, upper)
>> + colMeans(Xbig)} )
>>> describe(meansDiag[4,])
>> meansDiag[4, ]
>>   n  missing distinct Info Mean  Gmd  .05  .10  
>> .25
>> 1000  100135.23  0.0207435.2135.21
>> 35.22
>> .50  .75  .90  .95
>>   35.2335.2535.2635.26
>> 
>> lowest : 35.18360 35.19756 35.20098 35.20179 35.20622
>> highest: 35.26367 35.26635 35.26791 35.27251 35.27302

Re: [R] Generating samples from truncated multivariate Student-t distribution

[peter dalgaard]

It's not obvious to me that that marginal distribution of one component of a 
multivariate truncated t is the corresponding univariate truncated t.

In fact, I would expect it to differ because of tail-dependence effects, e.g.

> r <- rtmvt(1e5, c(30,0), diag(2), lower=c(29,-Inf), upper=c(31, +Inf), df=4)
> sd(r[,2])
[1] 1.191654
> r <- rtmvt(1e5, c(30,0), diag(2), lower=c(35,-Inf), upper=c(37, +Inf), df=4)
> sd(r[,2])
[1] 3.504233

-pd

> On 9 May 2017, at 19:09 , Czarek Kowalski  wrote:
> 
> Dear Members,
> I am working with 6-dimensional Student-t distribution with 4 degrees
> of freedom truncated to [20; 60]. I have generated 100 000 samples
> from truncated multivariate Student-t distribution using rtmvt
> function from package ‘tmvtnorm’. I have also calculated  mean vector
> using equation (3) from attached pdf. The problem is, that after
> summing all elements in one column of rtmvt result (and dividing by
> 100 000) I do not receive the same result as using (3) equation. Could
> You tell me, what is incorrect, why there is a difference?
> Yours faithfully
> Czarek Kowalski
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

-- 
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Office: A 4.23
Email: pd@cbs.dk  Priv: pda...@gmail.com

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Generating samples from truncated multivariate Student-t distribution

[Czarek Kowalski]

Previously I had used another language to make calculations based on
theory. I have calculated using R and I have received another results.
My theoretical calculation does not take into account the full
covariance matrix (only 6 elements from diagonal). Code based on
theory:

df = 4;   #degrees of freedom
sigmas = c(1, 4, 2, 5, 3, 6) # roots of diagonal elements of covariance matrix
meann = c(55, 40, 50, 35, 45, 30)
alfa1 = 20; # lower truncation
beta1 = 60; # upper truncation
a = (alfa1 - meann)/sigmas;
b = (beta1 - meann)/sigmas;
E = meann + sigmas * ((gamma(df - 1)/2)*((df + a^2)^(-(df-1)/2) - (df
+ b^2)^(-(df-1)/2))*df^(df/2))/(2*(pt(b,df)-pt(a,df))*gamma(df/2)*gamma(1/2))
E


Kind regards
Czarek

On 9 May 2017 at 23:50, David Winsemius  wrote:
>
>> On May 9, 2017, at 2:33 PM, David Winsemius  wrote:
>>
>>
>>> On May 9, 2017, at 2:05 PM, Czarek Kowalski  wrote:
>>>
>>> I have already posted that in attachement - pdf file.
>>
>> I see that now. I failed to scroll to the 3rd page.
>>
>>> I am posting
>>> plain text here:
>>>
 library(tmvtnorm)
 meann = c(55, 40, 50, 35, 45, 30)
 covv = matrix(c(  1, 1, 0, 2, -1, -1,
>>>   1, 16, -6, -6, -2, 12,
>>>   0, -6, 4, 2, -2, -5,
>>>   2, -6, 2, 25, 0, -17,
>>>  -1, -2, -2, 0, 9, -5,
>>>  -1, 12, -5, -17, -5, 36), 6, 6)
>>> df = 4
>>> lower = c(20, 20, 20, 20, 20, 20)
>>> upper = c(60, 60, 60, 60, 60, 60)
>>> X1 <- rtmvt(n=10, meann, covv, df, lower, upper)
>>>
>>>
 sum(X1[,1]) / 10
>>> [1] 54.98258
>>> sum(X1[,2]) / 10
>>> [1] 40.36153
>>> sum(X1[,3]) / 10
>>> [1] 49.83571
>>> sum(X1[,4]) / 10
>>> [1] 34.69571  # "4th element of mean vector"
>>> sum(X1[,5]) / 10
>>> [1] 44.81081
>>> sum(X1[,6]) / 10
>>> [1] 31.10834
>>>
>>> And corresponding results received using equation (3) from pdf file:
>>> [54.97,
>>> 40,
>>> 49.95,
>>> 35.31, #  "4th element of mean vector"
>>> 44.94,
>>> 31.32]
>>>
>>
>> I get similar results for the output from your code,
>>
>> My 100-fold run of your calculations were:
>>
>> meansBig <- replicate(100, {Xbig <- rtmvt(n=10, meann, covv, df, lower, 
>> upper)
>> colMeans(Xbig)} )
>>
>> describe(meansBig[4,])  # describe is from Hmisc package
>>
>> meansBig[4, ]
>>   n  missing distinct Info Mean  Gmd  .05  .10  
>> .25
>> 1000  1001 34.7  0.0195434.6834.68
>> 34.69
>> .50  .75  .90  .95
>>   34.7034.7234.7234.73
>>
>> lowest : 34.65222 34.66675 34.66703 34.66875 34.67566
>> highest: 34.72939 34.73012 34.73051 34.73742 34.74441
>>
>>
>> So agree, 35.31 is outside the plausible range of an RV formed with that 
>> package, but I don't have any code relating to your calculations from theory.
>
> Further investigation:
>
> covDiag <- covv*( row(covv)==col(covv) )  # just the diagonal means
>
> Repeat with all zero covariances:
>
>> meansDiag <- replicate(100, {Xbig <- rtmvt(n=10, meann, covDiag, df, 
>> lower, upper)
> + colMeans(Xbig)} )
>> describe(meansDiag[4,])
> meansDiag[4, ]
>n  missing distinct Info Mean  Gmd  .05  .10  
> .25
>  1000  100135.23  0.0207435.2135.21
> 35.22
>  .50  .75  .90  .95
>35.2335.2535.2635.26
>
> lowest : 35.18360 35.19756 35.20098 35.20179 35.20622
> highest: 35.26367 35.26635 35.26791 35.27251 35.27302
>
> So failing to account for the covariances in your theoretical calculations 
> mostly explains the apparent discrepancy, although your value of 35.31 would 
> be at the  far end of a statistical distribution and I wonder about some sort 
> of error in your theoretical calculation, which didn't appear to take into 
> account the covariance matrix.
>
> Best;
> David.
>
>
>
>>
>> Best;
>> David.
>>
>>
>>> On 9 May 2017 at 22:17, David Winsemius  wrote:

> On May 9, 2017, at 1:11 PM, Czarek Kowalski  
> wrote:
>
> Of course I have expected the difference between theory and a sample
> of realizations of RV's and result mean should still be a random
> variable. But, for example for 4th element of mean vector: 35.31 -
> 34.69571 = 0.61429. It is quite big difference, nieprawdaż? I have
> expected that the difference would be smaller because of law of large
> numbers (for 10mln samples the difference is quite similar).

 I for one have no idea what is meant by a "4th element of mean vector". So 
 I have now idea what to consider "big". I have found that my intuitions 
 about multivariate distributions, especially those where the covariate 
 structure is as complex as you have assumed, are often far from simulated 
 results.

 I suggest you post some code and results.

 --
 David.

Re: [R] Non-Linear Regression Help

[David Stevens]

I have a fair bit of experience with both nls and rating curves. This is 
not a nls() problem, this is a model problem. The power law rating curve 
favored by hydrologists would not apply to your data as it's based on 
the idea that a log-log plot of discharge vs. stage, or state+a in your 
case is a straight line, statistical assumptions notwithstanding. A 
log-log plot of your data,

plot(discharge~stage,data=yourdata,pch=19,log='xy')

clear is not a straight line. The very large discharge at stage=6.53 vs. 
stage=6.32 says one of two things: 1) there is an error in the data 
(perhaps the 2592.05 should be 592.05) or 2) the river channel geometry 
has changed dramatically, as in overtopping its banks or perhaps there's 
a smaller central channel set into a larger flood channel, similar to 
the LA river of the movies. The way forward is 1) recheck your data or 
2) recheck your data and use a two-piece model with one piece for stage 
<= 6.32 and a second piece for stage > 6.32. For this second approach to 
work, you'll need more data than you have given us here.

BTW, nls() should work fine if the model/data combination meet the 
requirements of 1) the model 'fits' the data, 2) the residuals are 
NIID(0,sigma^2), the parameters C, a, and n are identifiable from the 
data (should be if the last point is excluded). As always, you'll need 
good starting values for the parameters (get them from a log-log plot). 
You may find, based on the residuals, that linear regression (lm, glm) 
are most appropriate so that the errors meet the criteria of constant 
variance. If none of this makes sense, buy and study the book

Nonlinear regression analysis: Its applications, D. M. Bates and D. G. 
Watts, Wiley, New York, 1988. ISBN 0471-816434.

The nls() application is the easy part.


Good luck

David Stevens

On 5/4/2017 4:58 PM, Zachary Shadomy wrote:
> I am having some errors come up in my first section of code. I have no
> issue in plotting the points. Is there an easier method for creating a
> non-linear regression using C*(x+a)^n. The .txt file is named
> stage_discharge with the two variables being stage and discharge.
> The data is a relatively small file listed below:
>
> stage discharge
> 6.53 2592.05
> 6.32 559.5782
> 5.96 484.2151
> 4.99 494.7527
> 3.66 456.0778
> 0.51 291.13
>
>
>
>
>
>> power.nls<-nls(stage_dischargee$discharge~C*(stage_discharge$stage+a)^n,
> data=stage_discharge, start=list(C=4, a=0, n=1))
>> C<-coef(power.nls)["C"]
>> a<-coef(power.nls)["a"]
>> n<-coef(power.nls)["n"]
>> plot(stage_discharge$stage, stage_discharge$discharge, pch=17, cex=1.25,
> ylab='Discharge (cfs )', xlab='Stage (ft)', font.lab=2, main='Boone Creek\n
> St age-Discharge Curve')
>> curve(C*(x+a)^n, add=TRUE, col="red")
>   [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

-- 
David K Stevens, P.E., Ph.D.
Professor and Head, Environmental Engineering
Civil and Environmental Engineering
Utah Water Research Laboratory
8200 Old Main Hill
Logan, UT  84322-8200
435 797 3229 - voice
435 797 1363 - fax
david.stev...@usu.edu



[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Problem with choose.files(default=..., multi=FALSE)

[Keith Jewell]

Thanks for confirming that I wasn't being stupid :-}

When using default=pathlong I get the _correct_ starting directory...
(M:\test\Averyveryveryveryverylongfoldername\Averyveryveryveryverylongfoldername\Averyveryveryveryverylongfoldername) 

... both in the environment I indicated originally (Windows Server 2008 
R2 x64) and also in Windows 10 x64


Keith Jewell

On 09/05/2017 17:49, Duncan Murdoch wrote:

On 09/05/2017 12:06 PM, Keith Jewell wrote:

I'm very hesitant to suggest that there's a bug in such a venerable R
function, but I can't see what I'm doing wrong. Any comments are welcome


Yes, it looks like a bug.  One other thing I find a little strange: the
starting directory seems wrong when I have the pathlong default.  Did
you see that?  (I'm in Windows 10, not the same version as you.)

Duncan Murdoch



When using choose.files() where:
 default = something
 multi = FALSE
 selected file path is shorter than the default
... then the returned value is at least as long as the default,
characters from default appearing (wrongly) at the end of the returned
value.

Example, in which all but the first choose.files() select
"M:\\test\\target.dat". Note the last result.

 > pathlong <- choose.files(caption = "long")
 > pathlong # long file name to use as default for short selection
[1]
"M:\\test\\Averyveryveryveryverylongfoldername\\Averyveryveryveryverylongfoldername\\Averyveryveryveryverylongfoldername\\target.dat"

 > choose.files(caption = "short")  # no default without multi works
[1] "M:\\test\\target.dat"
 > choose.files(default=pathlong, caption = "short") # default without
multi= works
[1] "M:\\test\\target.dat"
 > choose.files(caption = "short", multi = FALSE) # multi = FALSE
without default works
[1] "M:\\test\\target.dat"
 > choose.files(default=pathlong, caption = "short", multi = TRUE) #
multi = TRUE with default works
[1] "M:\\test\\target.dat"
 > choose.files(default=pathlong, caption = "short", multi = FALSE) #
multi = FALSE with default fails
[1]
"M:\\test\\target.dat\\ryveryverylongfoldername\\Averyveryveryveryverylongfoldername\\Averyveryveryveryverylongfoldername\\target.dat"


 > # in case it's relevant
 > sessionInfo()
R version 3.4.0 (2017-04-21)
Platform: i386-w64-mingw32/i386 (32-bit)
Running under: Windows Server 2008 R2 x64 (build 7601) Service Pack 1

Matrix products: default

locale:
[1] LC_COLLATE=English_United Kingdom.1252  LC_CTYPE=English_United
Kingdom.1252
[3] LC_MONETARY=English_United Kingdom.1252 LC_NUMERIC=C

[5] LC_TIME=English_United Kingdom.1252

attached base packages:
[1] graphics  grDevices datasets  stats tcltk utils tools
   methods
[9] base

other attached packages:
  [1] CBRIutils_1.0   stringr_1.2.0   svSocket_0.9-57 TinnR_1.0-5
R2HTML_2.3.2
  [6] Hmisc_4.0-3 ggplot2_2.2.1   Formula_1.2-1   survival_2.41-3
lattice_0.20-35

loaded via a namespace (and not attached):
  [1] RColorBrewer_1.1-2  htmlTable_1.9   digest_0.6.12
htmltools_0.3.6
  [5] splines_3.4.0   scales_0.4.1grid_3.4.0
checkmate_1.8.2
  [9] devtools_1.12.0 knitr_1.15.1munsell_0.4.3
compiler_3.4.0
[13] tibble_1.3.0nnet_7.3-12 acepack_1.4.1
Matrix_1.2-10
[17] svMisc_0.9-70   plyr_1.8.4  base64enc_0.1-3
data.table_1.10.4
[21] stringi_1.1.5   magrittr_1.5gtable_0.2.0
colorspace_1.3-2
[25] foreign_0.8-68  cluster_2.0.6   gridExtra_2.2.1
htmlwidgets_0.8
[29] withr_1.0.2 lazyeval_0.2.0  backports_1.0.5
memoise_1.1.0
[33] rpart_4.1-11Rcpp_0.12.10latticeExtra_0.6-28
 >

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.





__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.