Maximum Likelihood Estimation for custom distribution

Hello,

I have the following problem.

I'm trying to estimate the parameters of a non-standard distribution. I launched dfittool, chose File -> Define Custom Distribution, created file with my own distribution and placed it to Documents/MATLAB/+prob/KouDistribution.m.

When I choose File -> Import Custom Distribution, I cann't find it in the list. I also copied my file to other different folders called '+prob'. But it had no result.

So, where should I place that file?

I have one more question.

I must write method fit for my new distribution, but I don't know, how. I now the probability density function, I know the cumulative density function.

May be, there is another way to estimate needed parameters?

This is a very important problem for me, so I ask someone to help me, who had the same problem, may be, or knows the solution.

Here is the code for distributions object

classdef KouDistribution < prob.ToolboxFittableParametricDistribution

    properties(Constant)
        DistributionName = 'kou';
    end

    properties(Dependent = true)
        mu
        sigma
        lambda
        p
        eta_1
        eta_2
    end

    properties(Constant)
        NumParameters = 6;
        ParameterNames = {'mu' 'sigma' 'lambda' 'p' 'eta_1' 'eta_2'};
        ParameterDescription = {'mu' 'sigma' 'lambda' 'p' 'eta_1' 'eta_2'};
    end

    properties(GetAccess = 'public', SetAccess = 'protected')
        ParameterValues
    end

    methods
        function pd = KouDistribution(mu, sigma, lambda, p, eta_1, eta_2)
            pd.ParameterValues = [mu sigma lambda p eta_1 eta_2];
            pd.ParameterIsFixed = [true true true true true true];
            pd.ParameterCovariance = zeros(pd.NumParameters);
        end

        function m = meat(this)
            m = this.mu;
        end
        function s = std(this)
            s = sqrt(2) * this.sigma;
        end
        function v = var(this)
            v = 2 * this.sigma^2;
        end
    end

    methods
        function this = set.mu(this, mu)
            this.ParameterValues(1) = mu;
        end
        function this = set.sigma(this, sigma)
            this.ParameterValues(2) = sigma;
        end
        function this = set.lambda(this, lambda)
            this.ParameterValues(3) = lambda;
        end
        function this = set.p(this, p)
            this.ParameterValues(4) = p;
        end
        function this = set.eta_1(this, eta_1)
            this.ParameterValues(5) = eta_1;
        end
        function this = set.eta_2(this, eta_2)
            this.ParameterValues(6) = eta_2;
        end

        function mu = get.mu(this)
            mu = this.ParameterValues(1);
        end
        function sigma = get.sigma(this)
            sigma = this.ParameterValues(2);
        end
        function lambda = get.lambda(this)
            lambda = this.ParameterValues(3);
        end
        function p = get.p(this)
            p = this.ParameterValues(4);
        end
        function eta_1 = get.eta_1(this)
            eta_1 = this.ParameterValues(5);
        end
        function eta_2 = get.eta_2(this)
            eta_2 = this.ParameterValues(6);
        end
    end

    methods(Static)
        function pd = fit(x, varargin)

        end
        function [nll, acov] = likefunc(params, x)
            mu = params(1);
            sigma = params(2);
            lambda = params(3);
            p = params(4);
            eta_1 = params(5);
            eta_2 = params(6);

            nll = -sum(log(pdffunc(x, mu, sigma, lambda, p, eta_1, eta_2)));
            acov = (sigma^2/n) * eye(2);
        end
        function y = cdffunc(x, mu, sigma, lambda, p, eta_1, eta_2)
            y = 0;
            for n = 1 : 10
                pi_n = exp(-lambda) * lambda^n / factorial(n);
                for k = 1 : n
                    pp = P(n, k, eta_1, eta_2, p);
                    qq = Q(n, k, eta_1, eta_2, p);
                    f1 = exp((sigma * eta_1)^2 / 2) / (sigma * sqrt(2 * pi)) ...
                        * (sigma * eta_1)^k * I(x - mu, -eta_1, -1 / sigma, -sigma * eta_1, k - 1);
                    f2 = exp((sigma * eta_2)^2 / 2) / (sigma * sqrt(2 * pi)) ...
                        * (sigma * eta_2)^k * I(x - mu, eta_2, 1 / sigma, -sigma * eta_2, k - 1);
                    y = y + pi_n * (pp * f1 + qq * f2);
                end
            end
            y = y + exp(-lambda) * normpdf(-(x - mu) / sigma);
        end
        function y = pdffunc(x, mu, sigma, lambda, p, eta_1, eta_2)
            y = 0;
            for n = 1 : 10
                pi_n = exp(-lambda) * lambda^n / factorial(n);
                for k = 1 : n
                    pp = P(n, k, eta_1, eta_2, p);
                    qq = Q(n, k, eta_1, eta_2, p);
                    h1 = Hh(-x / sigma + sigma * eta_1, k - 1);
                    h2 = Hh(x / sigma + sigma * eta_2, k - 1);
                    f1 = (sigma * eta_1)^k * exp((sigma * eta_1)^2 / 2) / (sigma * sqrt(2 * pi))...
                        * exp(-x * eta_1) .* h1;
                    f2 = (sigma * eta_2)^k * exp((sigma * eta_2)^2 / 2) / (sigma * sqrt(2 * pi))...
                        * exp(x * eta_2) .* h2;
                    y = y + pi_n * (pp * f1 + qq * f2);
                end
            end
            y = y + exp(-lambda) * normcdf(-(x - mu) / sigma);
        end
        function y = invfunc(p, mu, sigma, lambda, p, eta_1, eta_2)

        end
        function y = randfunc(mu, sigma, lambda, p, eta_1, eta_2, varargin)

        end
    end

    methods(Static, Hidden)
        function info = getInfo
            info.Name = 'Kou';
            info.code = 'kou';
            info.pnames = {'mu', 'sigma', 'lambda', 'p', 'eta_1', 'eta_2'};
            info.pdescription = {'mu', 'sigma', 'lambda', 'p', 'eta_1', 'eta_2'};
            info.prequired = [0, 0.2, 3, 0.4, 10, 5];
            info.hasconfbounds = false;
            info.censoring = false;
            info.support = [-Inf Inf];
            info.closedbound = [false false false false false false];
            info.iscontinuous = true;
            info.islocscal = true;
            info.uselogpp = false;
            info.optimopts = false;
            info.logci = [false false false false false false];
        end
    end
end

function f = Hh(x, n) f(1:length(x)) = 0; for i = 1:length(x) f(i) = integral(@(t) (t - x(i)).^n .* exp(-t.^2 / 2), x(i), Inf); end end

function f = P(n, k, eta_1, eta_2, p) q = 1 - p; f = 0; for i = k : n - 1 f = f + nchoosek(n - k - 1, i - k) * nchoosek(n, i) * (eta_1 / (eta_1 + eta_2))^(i - k)... * (eta_2 / (eta_1 + eta_2))^(n - i) * p^i * q^(n - i); end end

function f = Q(n, k, eta_1, eta_2, p) q = 1 - p; f = 0; for i = k : n - 1 f = f + nchoosek(n - k - 1, i - k) * nchoosek(n, i) * (eta_1 / (eta_1 + eta_2))^(n - i)... * (eta_2 / (eta_1 + eta_2))^(i - k) * p^(n - i) * q^(i); end end

function f = I(c, alpha, beta, delta, n) f(1:length(c)) = 0; for i = 1:length(c) f = integral(@(x) exp(alpha * x) .* Hh(beta * x - delta, n), c(i), Inf); end end

Answer 1

[Tom Lane]

Here are some tricks to help debug this sort of thing:

makedist -reset              % will cause MATLAB to revisit your file and re-load
p = makedist('kou',1,2,3...) % to make sure you can create one of these
pdf(p,1:4)                   % to make sure you can call its pdf method

Here are the things I had to change to get your file to load into dfittool:

1. The invfunc method had two inputs with the same name. Change one of them.

2. The getInfo method needs to start with

info = [e-mail]cDistribution('prob.KouDistribution');
info.name = 'Kou';

(You were missing the first line and had Name rather than name on the second.)

3. Correct the spelling of the islocscale property in this same method.

With these changes I could get dfittool to show your distribution. When I tried to fit it I got an error because your fit method isn't implemented yet.