3
$\begingroup$

I have a 720 hourly set of wind speed and wind direction data and I want to fit the Weibull Distribution on it. After searching for some time, I wrote the following code in Python to get my distribution, I will also share my image for clarification. windspeed weibull<code>enter code here</code>

ws.hist(bins=np.arange(0, ws.max()), alpha=0.5, normed=True)
weibull_params = (sc.stats.exponweib.fit(ws,floc=0,f0=1))
x= ws
def weib(x,lamb,k):
   return (k / lamb) * (x / lamb)**(k-1) * np.exp(-(x/lamb)**k)
k_shape, lamb_scale = weibull_params[1], weibull_params[3]

plt.plot(x, weib(x, lamb_scale, k_shape), label='self-defined Weibull')

Does this seem alright? How do I make sense of this? I also found a package here but I don't understand how to implement it to my code. https://github.com/cqcn1991/Wind-Speed-Analysis

Thank you.

$\endgroup$
  • $\begingroup$ Perhaps the data points are not in sequential order along the x-axis with is leading to a jump around of the line. $\endgroup$ – Edmund Jun 5 '19 at 21:38
  • $\begingroup$ It is actually a time series data. I used a histogram. Should I sort and then fit it? $\endgroup$ – cwanderroycbooks Jun 5 '19 at 21:55
  • $\begingroup$ I would check the order of the points in x along the x-axis as a first step. If they are not in sequential order along the x-axis then I would sort x and replot to see if it makes a difference. I suspect that it may. $\endgroup$ – Edmund Jun 5 '19 at 22:00
  • $\begingroup$ It stays the same actually. Maybe because I have too many bins. $\endgroup$ – cwanderroycbooks Jun 8 '19 at 5:06
2
$\begingroup$

Tried a few things in the OP comments that didn't work. However, you may apply Wolfram Language to your project. There is a free Wolfram Engine for developers and with the Wolfram Client Library for Python you can use these functions in Python.

import datetime

from wolframclient.evaluation import WolframLanguageSession
from wolframclient.language import wl, wlexpr

Start the Wolfram session

wolfSession = WolframLanguageSession()

I will collect my own WindSpeedData from a WeatherData weather station closest to the "Country" Entity of Bermuda during October 2016.

weatherStation = wolfSession.evaluate(
    wl.WeatherData(wl.Entity('Country','Bermuda'))
);
print(weatherStation)
Entity['WeatherStation', 'TXKF']
windData = wolfSession.evaluate(
    wl.WindSpeedData(
        weatherStation,
        [
            datetime.datetime(2016,10,1),
            datetime.datetime(2016,10,31)
        ]
    )
)

windData is a TimeSeries object with observation count comparable to your problem.

print(wolfSession.evaluate(windData('PathLength')))
1058

The result can be visualised with DateListPlot by Exporting in one of the supported Raster Image Formats or Vector Graphics Formats.

wolfSession.evaluate(
    wl.Export(
        '<path with image filename>', 
        wl.DateListPlot(windData, 
            PlotTheme='Detailed', 
            PlotRange=wl.All, 
            FrameLabel=wl.Automatic
        )
    )
)

Mathematica graphics

A hurricane passed over in that month with is the cause of that spike in wind speed.

Now that there is data Histogram on the PDF scale gives.

hist=wolfSession.evaluate(
    wl.Histogram(windData, wl.Automatic, 'PDF',
        PlotTheme='Detailed', 
        ChartStyle=wl.ColorData('Crayola','Silver'),
        PlotRange=wl.All
    )
);

wolfSession.evaluate(
    wl.Export(
        '<path with image filename>', 
        hist
    )
)

Mathematica graphics

Use EstimatedDistribution for the WeibullDistribution.

windDistrbution=wolfSession.evaluate(
    wl.EstimatedDistribution(
        wl.QuantityMagnitude(windData('Values')),
        wl.WeibullDistribution(wl.Global.alpha, wl.Global.beta, wl.Global.mu)
    )
);
print(windDistrbution) 
WeibullDistribution[1.883495945177254, 28.34295076324276, -0.7675654467340361]

Or use FindDistribution restricted to the continuous distribution functions to auto find and fit a distribution.

autoDistribution=wolfSession.evaluate(
    wl.FindDistribution(
        wl.QuantityMagnitude(windData('Values')),
        TargetFunctions='Continuous'
    )
);

print(autoDistribution)
ExtremeValueDistribution[18.436141153779957, 10.204118744677338]

Check for goodness of with DistributionFitTest.

print(
    wolfSession.evaluate(
        wl.Map(
            wl.Function(wl.DistributionFitTest(windData,wl.Slot(1),'PValue')), 
            [windDistrbution,autoDistribution]
        )
    )
)
[0.0002862317707539308, 0.1802410209978217]

The Weibull fit is not very good and Extreme Value not that strong. In any case the Plot the PDF functions.

pdfs=wolfSession.evaluate(
    wl.Map(
        wl.Function(wl.PDF(wl.Slot(1), wl.Global.x)), 
        [windDistrbution,autoDistribution]
    )
);


pdfPlot=wolfSession.evaluate(
    wl.Plot(
        pdfs,
        [wl.Global.x, 0, wl.QuantityMagnitude(wl.Max(windData))],
        PlotTheme='Detailed',
        PlotLegends=wl.Placed(
            wl.LineLegend(wl.Automatic, [windDistrbution,autoDistribution]), 
            wl.Below)
    )
)

wolfSession.evaluate(
    wl.Export(
        '<path with image filename>', 
        pdfPlot
    )
)

Mathematica graphics

Combine the plots with Show.

wolfSession.evaluate(
    wl.Export(
        '<path with image filename>', 
        wl.Show(hist, pdfPlot)
    )
)

Mathematica graphics

Terminate the session

wolfSession.terminate()

Hope this helps.

| improve this answer | |
$\endgroup$
  • $\begingroup$ Cool, didn‘t know that you can use Mathematica on Python... need to have a good look! $\endgroup$ – Peter Jun 8 '19 at 19:16
  • 1
    $\begingroup$ @Peter Its new since version 11.3. They are calling it "Wolfram Language" these days. It seems that Mathematica is slowing being rebranded as "Wolfram Desktop". $\endgroup$ – Edmund Jun 8 '19 at 19:20
  • 1
    $\begingroup$ Thanks, good to know. I worked with Mathematica some 7 years ago. I really liked it. Looking forward to try Wolfram Language! $\endgroup$ – Peter Jun 8 '19 at 19:39
  • $\begingroup$ Thanks a lot. I will try what you suggested. I have meanwhile attempted to fit a distribution to my dataset. $\endgroup$ – cwanderroycbooks Jun 18 '19 at 16:02
  • $\begingroup$ Thanks a lot. I will try what you suggested. I have meanwhile attempted to fit a distribution to my dataset. I actually get the following error wolframclient.exception.WolframKernelException: Cannot locate a kernel automatically. Please provide an explicit kernel path. $\endgroup$ – cwanderroycbooks Jun 18 '19 at 16:22

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.