Calculation of average damage using analytical methods

[someoneSomewhere]

TL;DR Summary

I have a task, the conditions of which you can read below. I understand how to solve it by approximating the value, but I would like to know how I can derive an equation (and is it even possible?), which can calculate the exact value.

Hello! A friend shared a problem he recently solved. It goes as follows:
Given:
Each dagger strike deals either normal damage = 20 or critical damage = 80.
After each strike, the probability of dealing normal and critical damage changes (initial probabilities are 90% and 10% respectively). The probability changes according to the following two conditions:
1. If the previous strike dealt normal damage, the probability of dealing normal damage decreases by 3%, and the probability of dealing critical damage increases by 3%.
2. If the previous strike dealt critical damage, the probability of dealing normal damage becomes 90%, and the probability of dealing critical damage becomes 10%.
Find:
The average damage for n strikes.

How to solve this problem through approximation is obvious to me. Here's the Python code:
The answer is approximately 31.5.

However, I am interested in the possibility of solving this problem analytically. That is, to write some equation that can be calculated for n values and get not an approximation, but an exact answer. Is this possible? I am weak in mathematical statistics and cannot imagine how to do this. You can simply point out some vectors of study that I need to explore to derive the equation myself, but I would be more grateful if you could derive this equation.

 

[willem2]

I can't derive an analytic expression. It is possible to write a recursive equation to compute the probability of critical damage after n rounds.
The process is a markov chain with 31 states.
State 0 is the starting state with critical hit probability 0.1. State 30 is the state with critical probability 1.
From state m you''ll go back to state 0 after a critical hit, and to state m+1 otherwise.
let prob(m,n) be the probability to be in state m after n rounds.
let crit(m) = 0.1 + 0.03m.
After 0 rounds, you'll be in the starting state so:
prob(0,0) = 1, and prob(m,0) = 0 if m <> 0;

prob(0,n) will be the sum of prob (m,n-1) * crit(m) (for all states m)
(you'll be in the starting state, if you got a critical hit in the round before)

prob(m,n) if m<>0 will be prob(m-1,n) * (1 - crit (m-1))
you'll be in state m, if you were in state m-1 on the round before and you dit not get a critical hit.

You can than sum prob(m, n) * crit(m) over all rounds to get the average damage.

from functools import cache

def crit(i):
    return 0.1 + 0.03 * i

@cache
def prob(m,n):

        if n == 0:
            if m == 0:
                return 1
            return 0

        if m == 0:
            probsum = 0
            for i in range (0,31):
                probsum += prob(i, n-1) * crit(i)
            return probsum

        return prob(m-1,n-1) * ( 1 - crit(m-1) )

crit_dam = 80
norm_dam = 20
rounds = 10

total_damage = 0
for n in range(0, rounds):
    for m in range(0,31):
        total_damage += prob(m, n)  * (crit(m) * crit_dam   +  (1-crit(m)) * norm_dam)

print (total_damage /rounds)

Your python program computes the average damage after 10000 stabs, and not after 10 stabs.
average damage after 1 stab is 26, after 10 30.103... and it seems to converge to about 31.471 if the rounds go to infinity.

This recursive program is horribly inefficient without using the @cache decorator from functools, and you'll probably never get an answer for n=20. If you want to implement it an another language you may have to make a 2 dimensional table for prob(m,n) yourself.

 

[someoneSomewhere]

I can't derive an analytic expression. It is possible to write a recursive equation to compute the probability of critical damage after n rounds.
The process is a markov chain with 31 states.
State 0 is the starting state with critical hit probability 0.1. State 30 is the state with critical probability 1.
From state m you''ll go back to state 0 after a critical hit, and to state m+1 otherwise.
let prob(m,n) be the probability to be in state m after n rounds.
let crit(m) = 0.1 + 0.03m.
After 0 rounds, you'll be in the starting state so:
prob(0,0) = 1, and prob(m,0) = 0 if m <> 0;

prob(0,n) will be the sum of prob (m,n-1) * crit(m) (for all states m)
(you'll be in the starting state, if you got a critical hit in the round before)

prob(m,n) if m<>0 will be prob(m-1,n) * (1 - crit (m-1))
you'll be in state m, if you were in state m-1 on the round before and you dit not get a critical hit.

You can than sum prob(m, n) * crit(m) over all rounds to get the average damage.

from functools import cache

def crit(i):
    return 0.1 + 0.03 * i

@cache
def prob(m,n):

        if n == 0:
            if m == 0:
                return 1
            return 0

        if m == 0:
            probsum = 0
            for i in range (0,31):
                probsum += prob(i, n-1) * crit(i)
            return probsum

        return prob(m-1,n-1) * ( 1 - crit(m-1) )

crit_dam = 80
norm_dam = 20
rounds = 10

total_damage = 0
for n in range(0, rounds):
    for m in range(0,31):
        total_damage += prob(m, n)  * (crit(m) * crit_dam   +  (1-crit(m)) * norm_dam)

print (total_damage /rounds)

Your python program computes the average damage after 10000 stabs, and not after 10 stabs.
average damage after 1 stab is 26, after 10 30.103... and it seems to converge to about 31.471 if the rounds go to infinity.

This recursive program is horribly inefficient without using the @cache decorator from functools, and you'll probably never get an answer for n=20. If you want to implement it an another language you may have to make a 2 dimensional table for prob(m,n) yourself.

Thank you very much for your answer. I think it more than satisfies my interest!

 

[pbuk]

I would like to know how I can derive an equation (and is it even possible?), which can calculate the exact value.

Yes you can calculate this analytically: the equation would involve some quite complicated summations which would be difficult to write out correctly. It is easier to write code to do it (see below), or to work it out using a spreadsheet.

You can simply point out some vectors of study that I need to explore to derive the equation myself

Just basic probability trees and a LOT of practice. Even with practice it is easy to make a mistake and I would always verify any result by Monte Carlo methods anyway.

How to solve this problem through approximation is obvious to me.

"Approximation" is not the right term here, we call the method you used "Monte Carlo simulation" or just "simulation".

Here is some code to calculate the answer and the results.

DAMAGE = 20
DAMAGE_CRIT = 80
P_CRIT_INITIAL = 0.1
P_CRIT_INCREMENT = 0.03

# Consecutive number of non-critical hits
k = 0

pLengths = []
pContinuing = 1
pCrit = P_CRIT_INITIAL
weightedSumOfLengths = 0

# Work around missing do...while in Python.
while True:
    # Calculate the probability of the sequence ending
    # here with a crit.
    pEnding = pContinuing * pCrit
    pContinuing -= pEnding
    pLengths.append(pEnding)

    # Accumulate the weighted sum of lengths.
    weightedSumOfLengths += (k + 1) * pEnding
    
    if pCrit >= 1: break
    
    # Mitigate truncation errors.
    pCrit = round(pCrit + P_CRIT_INCREMENT, 15)
    k += 1

nonCrits = 0
weightedSum = 0
totalProbability = 0

formatter = ' {:3} {:10.4g} {:10.6f}'
print('Hits Probability   Damage')

# Iterate over all possible sequence lengths.
for pLength in pLengths:
    # The probability that a hit is within a sequence of length nonCrits + 1.
    pInSequenceOfLength = pLength * (nonCrits + 1) / weightedSumOfLengths
    # Sum the probabilities for checking later.
    totalProbability += pInSequenceOfLength
    damagePerHit = (nonCrits * DAMAGE + DAMAGE_CRIT) / (nonCrits + 1)
    # Aggregate the weighted damage
    weightedSum += damagePerHit * pInSequenceOfLength

    nonCrits += 1
    print(formatter.format(nonCrits, pInSequenceOfLength, weightedSum))

print('Total probability (should be 1)', sum(pLengths), totalProbability)
print(weightedSum)

Hits Probability   Damage
   1    0.01912   1.529585
   2    0.04474   3.766604
   3    0.07186   6.641000
   4    0.09557   9.986079
   5      0.112  13.571601
...
Total probability (should be 1) 0.9999999999999999 0.9999999999999998
31.471888847439825