[R] Need help..

[Ajit K Jena]

Dear All,

I am trying to generate some Pareto random variates
using the inverse method. This is really straightfoward
and my R function looks as below :

pareto <- function(c, a, cnt=1000)
{
	u <- runif(cnt)

	x <- (c / ((u ^ (1 / a))))

	mean.theo <- ((c * a) / (a - 1))
	mean.gen <- mean(x)
	cat('Pareto mean : theoritical', mean.theo, 
		'generated', mean.gen, '\n') 

	return(x)
}

I am also giving the outputs below to illustrate my
point. I understand the point that as 'a' is less
than 2.0, the variance does not exist. But as 'a'
approaches 1.0, the mean goes down approximately to
1/3 of the theoritical mean. I understand it is
getting slowly into unstable mean syndrome but I am
really baffled that it should remain so low consistently.

Am I missing something here ? Can anyone please suggest
some better method ? 

Regards.

--ajit

> pareto(1000, 1.6, 5000) -> k
Pareto mean : theoritical 2666.667 generated 2533.681 
> pareto(1000, 1.6, 5000) -> k
Pareto mean : theoritical 2666.667 generated 3182.131 
> pareto(1000, 1.6, 5000) -> k
Pareto mean : theoritical 2666.667 generated 2696.346 
> pareto(1000, 1.6, 5000) -> k
Pareto mean : theoritical 2666.667 generated 2624.048 
> pareto(1000, 1.6, 5000) -> k
Pareto mean : theoritical 2666.667 generated 2465.958 
> pareto(1000, 1.6, 5000) -> k
Pareto mean : theoritical 2666.667 generated 2968.517 

There are some big swings in the set above, but still ok I think... 

> pareto(1000, 1.4, 5000) -> k
Pareto mean : theoritical 3500 generated 3162.978 
> pareto(1000, 1.4, 5000) -> k
Pareto mean : theoritical 3500 generated 3129.814 
> pareto(1000, 1.4, 5000) -> k
Pareto mean : theoritical 3500 generated 3203.407 
> pareto(1000, 1.4, 5000) -> k
Pareto mean : theoritical 3500 generated 3379.751 
> pareto(1000, 1.4, 5000) -> k
Pareto mean : theoritical 3500 generated 3345.743 
> pareto(1000, 1.4, 5000) -> k
Pareto mean : theoritical 3500 generated 3387.779 

The above are still sensible I would think...

> pareto(1000, 1.2, 5000) -> k
Pareto mean : theoritical 6000 generated 4335.049 
> pareto(1000, 1.2, 5000) -> k
Pareto mean : theoritical 6000 generated 4929.394 
> pareto(1000, 1.2, 5000) -> k
Pareto mean : theoritical 6000 generated 4806.154 
> pareto(1000, 1.2, 5000) -> k
Pareto mean : theoritical 6000 generated 7470.355 
> pareto(1000, 1.2, 5000) -> k
Pareto mean : theoritical 6000 generated 5078.184 
> pareto(1000, 1.2, 5000) -> k
Pareto mean : theoritical 6000 generated 4861.391 
> pareto(1000, 1.2, 5000) -> k
Pareto mean : theoritical 6000 generated 11153.21 
> pareto(1000, 1.2, 5000) -> k
Pareto mean : theoritical 6000 generated 5626.684 

The trend has started above ..

> pareto(1000, 1.1, 5000) -> k
Pareto mean : theoritical 11000 generated 5752.399 
> pareto(1000, 1.1, 5000) -> k
Pareto mean : theoritical 11000 generated 7003.196 
> pareto(1000, 1.1, 5000) -> k
Pareto mean : theoritical 11000 generated 6012.118 
> pareto(1000, 1.1, 5000) -> k
Pareto mean : theoritical 11000 generated 5761.349 

Things below are really nasty...

> pareto(1000, 1.06, 5000) -> k
Pareto mean : theoritical 17666.67 generated 5897.83 
> pareto(1000, 1.06, 5000) -> k
Pareto mean : theoritical 17666.67 generated 5833.972 
> pareto(1000, 1.06, 5000) -> k
Pareto mean : theoritical 17666.67 generated 6125.007 
> pareto(1000, 1.06, 5000) -> k
Pareto mean : theoritical 17666.67 generated 7844.268 
> pareto(1000, 1.06, 5000) -> k
Pareto mean : theoritical 17666.67 generated 6741.457 
> pareto(1000, 1.06, 5000) -> k
Pareto mean : theoritical 17666.67 generated 9166.882 
> pareto(1000, 1.06, 5000) -> k
Pareto mean : theoritical 17666.67 generated 6726.428 

--------------------------------------------------------------------------
Ajit K. Jena					  Phone  :  +46-455-385655
						  Fax    :  +46-455-385667
Dept. of Telecom and Signal Processing		  Mobile :  +46-708-876803
University of Karlskrona/Ronneby		  Email  : akj at its.hk-r.se
S-371 79 Karlskrona, Sweden
--------------------------------------------------------------------------

On Tue, 23 Nov 1999, T. V. Kurien wrote:

> Dear akj,
>     The inversion method IS THE ONLY WAY (in general) to generate
> rv from a particular CDF. For one thing, the CDF := P(X >=x) = F(x). For
> the pareto, it is as you have stated, except, of course that x>c. In terms
> of the mean, it is ac/(a-1) for a>1 .. I assume that you are aware of this,
> and on the strong dependence on c.
> 
> Let me re-iterate: The inverse method is CERTAINLY GOOD ENOUGH. Make sure
> that you have not screwed up in the random variate generation (inversion
> itself), in this case, F(x) is pretty trivial to invert.
> v
> 


-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._