Calculator of skill-dependent power and marginal contribution of collective dichotomous choice (Ben-Yashar, Nitzan, and Tajika, 2021; Nitzan and Tajika, 2021).

using LinearAlgebra, Statistics

using Random

#Calculation for Marginal Contribution

#Calculate the optimal-weight profile for given skill profile

function OpW(skills)

    weights = similar(skills)

    for (i,s) in enumerate(skills)

        weights[i]= log(s/(1-s))

    end

    return weights

end    

#Decision  patterns

function Dec(dim)

    Dec = zeros(2^dim, dim)

    for i in 1:dim

        CC = vcat(-ones(2^(dim-i)),ones(2^(dim-i)))

        Dec[:, i]= repeat(CC, 2^(i-1))

    end

    return Dec

end

#Positive indicator

function cordec(s)

    if s > 0

        return 1

    else

        return 0

    end

end

#b or 1-b

function binarypr(a,b)

    if a == 1

        return b

    else

        return 1-b

    end

end

#Calculating the probability of correct decision making for given skill vector

function prob_correct_decision(skills)

    dim = size(skills,2)

    weight = OpW(skills)

    decision = zeros(2^dim)

    prob = zeros(2^dim)

    indv = Dec(dim)

    for i in 1:2^dim

        decision[i] = cordec(dot(indv[i,:],weight))

        prob[i] = prod(binarypr(indv[i,j],skills[j]) for j in 1:dim)

        i += 1

    end

    return dot(decision,prob)

end

#calculating marginal contribution

function marginal_contribution(skills)

    dim = size(skills,2)

    maxpr = prob_correct_decision(skills)

    mc_v = zeros(dim)

    for i in 1:dim

        others = deleteat!(skills[1,:],i)'

        mc_v[i] = maxpr - prob_correct_decision(others)

        i += 1

    end

    return mc_v

end

#Calculating s-d power

function sd_power(skills)

    dim = size(skills,2)

    weight = OpW(skills)

    decision = zeros(2^dim)

    ndecision = zeros(2^dim)

    prob = zeros(2^dim)

    indv = Dec(dim)

    sd = zeros(dim)

# k=1    

        nweight=copy(weight)

        nweight[1]=0

        for i in 1:2^dim

            decision[i] = cordec(dot(indv[i,:],weight))

            ndecision[i] = cordec(dot(indv[i,:],nweight))

            prob[i] = prod(binarypr(indv[i,j],skills[j]) for j in 2:dim)

            sd[1] = sd[1]+ prob[i]*abs(decision[i] - ndecision[i] )

            i += 1

        end

# k = 2:dim-1    

    for k in 2:dim-1

        nweight=copy(weight)

        nweight[k]=0

        for i in 1:2^dim

            decision[i] = cordec(dot(indv[i,:],weight))

            ndecision[i] = cordec(dot(indv[i,:],nweight))

            prob[i] = prod(binarypr(indv[i,j],skills[j]) for j in 1:k-1)*prod(binarypr(indv[i,j],skills[j]) for j in k+1:dim)

            sd[k] = sd[k]+ prob[i]*abs(decision[i] - ndecision[i] )

            i += 1

        end

    k += 1

    end

# k=dim    

        nweight=copy(weight)

        nweight[dim]=0

        for i in 1:2^dim

            decision[i] = cordec(dot(indv[i,:],weight))

            ndecision[i] = cordec(dot(indv[i,:],nweight))

            prob[i] = prod(binarypr(indv[i,j],skills[j]) for j in 1:dim-1)

            sd[dim] = sd[dim]+ prob[i]*abs(decision[i] - ndecision[i] )

            i += 1

        end    

    return sd

end

#Gini index

function mygini(dist)

    n=size(dist,1)

    gini = (1- sum(sum(minimum((dist[i],dist[j])) for j in 1:n) for i in 1:n)/((n^2)*mean(dist)))

    return gini

end

#Simulation with gini index

#N: Number of members, T: Number of trials.

function simulation_gini(N,T) 

    gini_m=zeros(T)

    gini_s = zeros(T)

    gini_w = zeros(T)

    for i in 1:T

        skills = rand(N)'./2 .+0.5

        gini_m[i] = mygini(marginal_contribution(skills))

        gini_s[i] = mygini(sd_power(skills))

        gini_w[i] = mygini(OpW(skills'))

        i += 1

    end

    return mean(gini_m), mean(gini_s), mean(gini_w)

end

# Coefficient of variation 

function cv(x)

    if mean(x) == 0

        return 0

        else 

    return std(x,corrected=false)/(mean(x)*sqrt(size(x)[1]-1))

    end

end

function simulation_cv(N,T)

    gini_m=zeros(T)

    gini_s = zeros(T)

    gini_skill = zeros(T)

    for i in 1:T

        skills = rand(N)'./2 .+0.5

        gini_m[i] = cv(marginal_contribution(skills))

        gini_s[i] = cv(sd_power(skills))

        gini_skill[i] = cv(OpW(skills'))

        i += 1

    end

    return mean(gini_m), mean(gini_s), mean(gini_skill)

end