qpowerlaw[p_, alpha_, xmin_, lower_, logp_] := {
    lower = True;
    logp = False;
    If[lower == False, p = -1 - p];
    If[logp, -Exp[p]];
    xmin*(1 - p)^(-1/(alpha - 1));

rpowerlaw[n_, xmin_, alpha_] := { 
    powList = Table[qpowerlaw[RandomVariate[UniformDistribution[{0,1}]]], {i,0, n}]; 

I am trying to use these two functions to properly write these lines of code in R, but I can't figure out how quite to do it.

Here is the code in R:

qpowerlaw <- function(p, alpha=2, xmin=1, lower.tail=T, log.p = F) {
     if (!lower.tail)
         p <- 1-p
     if (log.p)
         p <- exp(p)
     xmin * ((1 - p) ^ (-1 / (alpha - 1)))
rpowerlaw <- function(n, alpha=2, xmin=1) {
     qpowerlaw(runif(n, 0, 1), alpha, xmin)


I want to be able to test alpha = 1.8 and xmin = 10, but when I run my code. I am getting results in powlist that are less than one, which can not be right. I tried to use ParetoDistribution to properly create the data, but I could not produce the kind of results i wanted (the alpha values between the ParetoDistribution and PowerLaw are very different). I also tried PowerDistribution, which lacks the xmin parameter, so I couldn't use that function either. If anyone knows, something I could try or any errors in my code. It would be very much appreciated.

  • $\begingroup$ How are you calling your function? One obvious problem I see is that qpowerlaw[] requires 5 variables, but in rpowerlaw[] you are calling it with just one. $\endgroup$
    – Feyre
    Dec 5 '16 at 15:28
  • $\begingroup$ Thank you, I didn't notice that, and it fixed up the issue. $\endgroup$ Dec 5 '16 at 20:29

You need to check the typos in your code. (Sometimes you are translating "<-" in R to "=-" in Mathematica when it should translate just to "=".) Also, the R code results in a vector of length n rather than n+1 which is what you get from {i,0,n}.

Does the following do what you want?

qpowerlaw[p_, alpha_: 2, xmin_: 1, lower_: True, logp_: False] := {
  If[! lower, p = 1 - p];
  If[logp, p = Exp[p]];
  xmin*(1 - p)^(-1/(alpha - 1))}

rpowerlaw[n_, alpha_: 2, xmin_: 1] :=
  Table[qpowerlaw[RandomVariate[UniformDistribution[{0, 1}]], alpha, xmin][[1]], {i, n}]

rpowerlaw[12, 2, 10]
(* {322.81124530346244`,48.06374170498156`,18.70766795287194`,18.129813319193655`,
10.103780695831814`,24.197428848689704`,44.91040965366457`,18.41061541949269`} *)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.