On Wed, Feb 01, 2012 at 11:01:21AM -0500, Val wrote:
> Hi Petr,
>
> Thank you very much. It looks we are almost there.
>
>
> #an example, use a function approximating your graphs
> * f <- function(x) { 0.7 - 1.25*((x - 1)^2 - 0.4)^2 }
> f <- function(x) { 0.2 - 1.5*((x - 1)^2 - 0.1)^2 }
>
> #plot function f(x)
> # x <- seq(0, 2, length=51)
> x <- seq(0.01, 1.75, length=51)
>
> y <- f(x)
> plot(x, y, type="l")
>
> The above plot may give the desired result.
>
> Here are the the conditions.
>
> 1. Y- value is always positive.
I used a polynomial as f(x), since this is simple.
However, for a better tuning of the shape, a spline
may be better.
> 2. Can I set the mean and SD of the value of X?
> Example mean of x= 24 and SD =5
For these conditions, it may be better to shift the
function by 24 to the left, so the mean will be 0
and the second moment will be 25. Then try to find
a spline with a shape, which you require, and with
first moment 0 and second moment 25 computed as follows.
# example spline
library(splines)
xOrig <- 2*(-5:5)
yOrig <- c(0, 3, 3.9, 4, 3.9, 3.9, 3.9, 4, 3.9, 3, 0)
ispl <- interpSpline(xOrig, yOrig)
# see the shape
x <- seq(-10, 10, length=10001)
y <- predict(ispl, x)$y
stopifnot(y >= 0)
plot(x, y, type="l")
# compute the moments
moment1 <- sum(x*y)/sum(y)
moment2 <- sum(x^2*y)/sum(y)
Here, we get
moment1
[1] -2.484611e-17
moment2
[1] 25.88646
I think that a trial and error method may be sufficient
for adjusting the original points to get [0, 25]
within a small error.
Another option is to create several candidate splines
and consider their mixture. The moments are linear
functions of the mixture parameters, so the mixture
parameters may be obtained as a solution of a system
of linear equations. If we have three splines such
that the triangle between the corresponding points
[moment1, moment2] contains the point [0, 25], then
the mixture parameters will be nonnegative and the mixture
well defined.
Petr.
______________________________________________
R-help@r-project.org mailing list
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.
