Project-Platypus / Platypus

A Free and Open Source Python Library for Multiobjective Optimization
GNU General Public License v3.0
547 stars 151 forks source link

Platypus for creating histograms with known mean, CV, skewness, and kurtosis #93

Closed tgmueller closed 1 year ago

tgmueller commented 5 years ago

Hi Folks:

And I thought a platypus was a beaver with duck feet..... Thanks @dhadka for the hard work you put into developing this.

Actually, I wrote code to try to obtain a distribution with a target mean, variance, skewness, and kurtosis. I almost have it working but I am having some problems. I'm using scipy.stats.gengamma to generate the distributions. I'm trying different weighting to adjust for the variable scale. It seems like I should have four constraints: problem = {i.e., Problem(4, 4, 4)} but I get an error. I am able to set the problem.constraints to "==0". If anyone has time to chat about this, please contact me at 515 745 7588 or tgmueller1@gmail.com.

from IPython.html.widgets import *

# from __future__ import print_function
from ipywidgets import IntSlider, link, VBox
from ipywidgets import interact, interactive, fixed, interact_manual
from ipywidgets import *
import ipywidgets as widgets
import numpy as np
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
import ipywidgets as widgets
import scipy.optimize
import scipy.stats
import csv
import pprint
import math
from platypus import *
import matplotlib.pyplot as plt
from scipy import interpolate
import numpy as np
from math import pow

def GammaDistribution(x):
    a = x[0]
    c = x[1]
    loc = x[2]
    scale = x[3]
    distribution = scipy.stats.gengamma(a=a, c=c, loc=loc, scale=scale)
    type(distribution)
    sample = distribution.rvs(size=4000)
    average_ = np.mean(sample)
    variance_ = np.var(sample)
    toup = distribution.stats("mvsk")
    average = str(np.asscalar(toup[0]))
    variance = str(np.asscalar(toup[1]))
    skewness = str(np.asscalar(toup[2]))
    kurtosis = str(np.asscalar(toup[3]))
    cv = math.sqrt(float(variance)) / float(average) * 100
    list_ = [
        (
            a,
            c,
            loc,
            scale,
            float(average),
            float(variance),
            cv,
            float(skewness),
            float(kurtosis),
        )
    ]
    listing.append((list_))
    #     print(
    #         "mean = ",
    #         str(average),
    #         ", variance = ",
    #         str(variance),
    #         "CV = ",
    #         str(cv),
    #         ", skewness = ",
    #         str(skewness),
    #         ", kurtosis = ",
    #         str(kurtosis),
    #     )
    #     print(
    #         "a = ", str(a), ", c = ", str(c), "loc = ", str(loc), ", scale = ", str(scale)
    #     )
    Dm = pow((Tm - float(average)), 5)
    Dcv = pow(Tcv - float(cv), 2)
    Ds = pow((Tk - float(skewness)) * 5, 3.0)
    Dk = (Ts - float(kurtosis)) * 10
    #     MinimizeThis = Dm + Dcv + Dk + Ds
    #     print(MinimizeThis)
    #     MinimizeThis = int(MinimizeThis * 10)
    #     print(MinimizeThis)
    return Dm, Dcv, Ds, Dk

#     return float(average), float(cv), float(skewness), float(kurtosis)
Tm = 22
Tcv = 10
Ts = -0.43
Tk = 0.352
# [12.328607937460614, -25.460218010291722, -199.09204975419763, 245.10235332412498] [-4.5894780344442605, 0.7954700620012862, 1.68761441185389e-05, -0.6432981836591487]
# [3.0351398891880197, -34.01527671333063, -108.83656356317944, 133.5798973594887] [0.4772643635389579, 1.8190410981758542, -0.03635086594101655, -1.3640100400096238]
# [5.271454107728743, -18.36716651547665, -90.48574624661777, 121.52803365906294] [0.4330459883351526, 10.28775634075794, -0.006536379702673841, -1.0141896837243147]
# [96.32804231700625, -8.22562709871211, -96.75381176324933, 207.24023880854904] [-0.018560531120510124, 3.356884702282528, 0.0450856297572932, -0.4687654547552718]
# [33.42107444373721, -25.203226221581353, -155.47004882428635, 203.62193078506584] [0.04174970934747253, 4.372548788239408, 0.15664124538884075, -0.5062796970829704]
# problem = Problem(4, 4, 4)
problem = Problem(4, 4)
problem.types[:] = [Real(1, 100), Real(-100, 50), Real(-200, 0), Real(0.2, 300)]
problem.constraints[:] = "==0"
problem.function = GammaDistribution
algorithms = [
    NSGAII,
    (NSGAIII, {"divisions_outer": 12}),
    (CMAES, {"epsilons": [0.05]}),
    GDE3,
    IBEA,
    (MOEAD, {"weight_generator": normal_boundary_weights, "divisions_outer": 12}),
    (OMOPSO, {"epsilons": [0.05]}),
    SMPSO,
    SPEA2,
    (EpsMOEA, {"epsilons": [0.05]}),
]
approach = GDE3
algorithm = approach(problem)
algorithm.run(200)
# print(algorithm.result)
print(Tm, Tcv, Ts, Tk)

for solution in unique(nondominated(algorithm.result)):
    print(solution.variables, solution.objectives)
    solution.objectives[0]
    print(test)
tgmueller commented 5 years ago

I also posted this question on stackoverflow: https://stackoverflow.com/questions/55541238/platypus-for-creating-histograms-with-known-mean-variance-skewness-and-kurtos

github-actions[bot] commented 1 year ago

This issue is stale and will be closed soon. If you feel this issue is still relevant, please comment to keep it active. Please also consider working on a fix and submitting a PR.