text: #!/usr/bin/env python3
  # -*- coding: utf-8 -*-
  """
  Created on Wed Feb 17 12:20:49 2021
  
  @author: atoosa
  """
  import numpy as np
  
  # uniform distribution
  def uniformDistr(a, b, size):
      data = np.random.uniform(low=a, high=b, size=size)
      return data
  
  def uniformMoments(a, b):
      m1 = (a+b)/2
      m2 = ((b-a)**2)/12
      return m1, m2
  
  # exponential distribution
  def expDistr(beta, size):
      data = np.random.exponential(scale=beta, size=size)
      return data
  
  def expMoments(beta):
      m1 = beta
      m2 = beta ** 2
      return m1, m2
  
  # normal distribution
  def normalDistr(mue, sigma, size):
      data = np.random.normal(loc=mue, scale=sigma, size=size)
      return data
  
  def normalMoments(mue, sigma):
      m1 = mue
      m2 = sigma ** 2
      return m1, m2
  
  # power-law distribution
  # https://stats.stackexchange.com/questions/173242/random-sample-from-power-law-distribution
  # https://arxiv.org/pdf/0706.1062.pdf : equation D4
  # using inverse of CDF
  def powerlawDistr(gamma, a, size):
      r = np.random.uniform(low=0, high=1, size=size)
      data = a * (1-r)**(-1/(gamma-1))
      return  data
  
  def powerlawMoments(gamma, a):
      if gamma>2:
          m1 = a * ((1-gamma)/(2-gamma))
      else:
          m1 = np.inf
      if gamma>3:
          m2 = a**2 * (gamma-1) * ((1/(gamma-3))-((gamma-1)/((gamma-2)**2)))
      else:
          m2 = np.inf
      
      return m1, m2
  
  np.random.seed(10)
  
  # uniform
  # calculate sample moments
  print("uniform:")
  a = 1
  b = 10
  for N in [1, 10, 100, 1000, 10000]:
      print("N: " + str(N))
      sample_m1 = []
      sample_m2 = []
      for K in range(0, 30):
          data = uniformDistr(a, b, N)
          sample_m1.append(np.mean(data))
          sample_m2.append(np.var(data))
      m1 = np.mean(sample_m1)
      m2 = np.mean(sample_m2)
      print("sample moments: "+ str(m1) + " ,"+ str(m2))
      
  # calculate moments from formula    
  true_m1, true_m2  = uniformMoments(a, b)
  print("calculated moments: "+ str(true_m1) + " ,"+ str(true_m2))
  
  # expoonential
  # calculate sample moments
  print("exponential:")
  beta = 1
  for N in [1, 10, 100, 1000, 10000]:
      print("N: " + str(N))
      sample_m1 = []
      sample_m2 = []
      for K in range(0, 30):
          data = expDistr(beta, N)
          sample_m1.append(np.mean(data))
          sample_m2.append(np.var(data))
      m1 = np.mean(sample_m1)
      m2 = np.mean(sample_m2)
      print("sample moments: "+ str(m1) + " ,"+ str(m2))
      
  # calculate moments from formula    
  true_m1, true_m2  = expMoments(beta)
  print("calculated moments: "+ str(true_m1) + " ,"+ str(true_m2))
  
  # normal
  # calculate sample moments
  print("normal:")
  mue = 1
  sigma = 10
  for N in [1, 10, 100, 1000, 10000]:
      print("N: " + str(N))
      sample_m1 = []
      sample_m2 = []
      for K in range(0, 30):
          data = normalDistr(mue, sigma, N)
          sample_m1.append(np.mean(data))
          sample_m2.append(np.var(data))
      m1 = np.mean(sample_m1)
      m2 = np.mean(sample_m2)
      print("sample moments: "+ str(m1) + " ,"+ str(m2))
      
  # calculate moments from formula    
  true_m1, true_m2  = normalMoments(mue, sigma)
  print("calculated moments: "+ str(true_m1) + " ,"+ str(true_m2))
  
  # powerlaw
  # calculate sample moments
  print("power-law:")
  a = 1
  gamma = 4
  for N in [1, 10, 100, 1000, 10000]:
      print("N: " + str(N))
      sample_m1 = []
      sample_m2 = []
      for K in range(0, 30):
          data = powerlawDistr(gamma, a, N)
          sample_m1.append(np.mean(data))
          sample_m2.append(np.var(data))
      m1 = np.mean(sample_m1)
      m2 = np.mean(sample_m2)
      print("sample moments: "+ str(m1) + " ,"+ str(m2))
      
  # calculate moments from formula    
  true_m1, true_m2  = powerlawMoments(gamma, a)
  print("calculated moments: "+ str(true_m1) + " ,"+ str(true_m2))
  