thr3ads.net - R help - [R] Need help.. [Nov 1999]

If this information is useful, please help other people find it:
Share via:

Ajit K Jena

1999-Nov-24 17:14 UTC

[R] Need help..

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) -> kPareto 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) -> kPareto 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) -> kPareto 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) -> kPareto 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) -> kPareto 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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

Prof Brian D Ripley

1999-Nov-24 19:36 UTC

head link

Means of pareto samples ([R] Need help..)

On Wed, 24 Nov 1999, Ajit K Jena wrote:
> 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 ? 
You need to be more patient. The mean is dominated by a few 
very large and very rare values, and those values are always 
positive.  The 99.99% point of your example with a=1.1 is about
4.5e6. Consider
> pareto(1000, 1.1, 50000) -> junk[1] 141980763
Pareto mean : theoritical 11000 generated 10499.47 > pareto(1000, 1.1, 50000) -> junk[1] 19794974
Pareto mean : theoritical 11000 generated 6782.521 

where I am printing the maximum values generated. The sample mean is just
extremely variable, with a very long right tail to its distribution too.

What I really mean is that you can't expect to estimate the mean
by simulation. Hopefully what you really want this for is less
sensitive. It is possible that you will see problems with poor
PRNGs in studying extremal behaviour, but I no evidence that is the 
issue here.


-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272860 (secr)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

Maybe Matching Threads

Search for more maybe matching threads

R help - Nov 1999 - Need help..

[R] Need help..

Means of pareto samples ([R] Need help..)

Maybe Matching Threads